1 Data wrangling

Load raw data

#stock assessment data shared by NMFS
#northeast region
ne_data <- readxl::read_excel(file.path(dir_anx, '_raw_data/NOAA_NMFS/stock_assessments/NE_Assessment_Records.xlsx'))
#midatlantic region
midatl_data <- readxl::read_excel(file.path(dir_anx, '_raw_data/NOAA_NMFS/stock_assessments/MA_Assessment_Records.xlsx'))

Clean data

#select just the columns we are interested in

data <- ne_data %>%
  bind_rows(midatl_data) %>%
  select(year = Year, 
         stock = Stock, 
         f_fmsy = `F/Fmsy`, 
         b_bmsy = `B/Bmsy`) %>%
  gather(key = "indicator", value = "value", -year, -stock) %>%
  filter(value >= 0) %>% #Smooth skate has F/Fmsy value of -9999 so removing that
  mutate(stock = str_replace(stock, " \\| Asmt & Status", ""),
         stock = str_replace(stock, " \\| Asmt", ""))
  
ggplot(data, aes(x = year, y = value, color = indicator)) +
  geom_hline(yintercept = 1, color = "darkgray") +
  geom_line() +
  facet_wrap(~stock, labeller = labeller(stock = label_wrap_gen(width = 25)), scales = "free_y") +
  theme_bw() +
  theme(strip.text = element_text(size=8))

Save data

write.csv(data, file = "data/nmfs_stock_assessment_data.csv")

Span plot

stock_ids <- data$stock %>% unique()

spanplot_df <- data %>%
  group_by(stock) %>%
  mutate(yr_min = min(year), yr_max = max(year))

span_plot <- ggplot(spanplot_df, aes(x = stock)) +
  #ggtheme_basic +
  theme(panel.border     = element_blank(),
        panel.grid.minor = element_blank(),
        panel.grid.major = element_line(colour = 'grey90'),
        panel.background = element_blank()) +
  geom_linerange(aes(ymin = yr_min, ymax = yr_max, 
                     color = indicator),
                 position = position_dodge(width = .5), alpha = .8) +
  scale_color_manual(values = c('darkred', 'darkblue'),
                       labels=c("B/Bmsy", "F/Fmsy")) +
  coord_flip() +
  labs(x = 'Stock',
       y = 'Year',
       color = 'Indicator')

span_plot

What species have multiple stocks? These ones will give us a tough time when matching up catch to the stocks.

dups <- data %>%
  separate(stock, into = c("species", "location"), sep = " - ") %>%
  select(-indicator, -value) %>%
  distinct() %>%
  group_by(year, species) %>%
  mutate(count = n()) %>%
  ungroup() %>%
  filter(count > 1) %>%
  mutate(stock = paste0(species, " - ", location)) %>%
  select(stock, species, location) %>%
  distinct()

dups

## # A tibble: 20 x 3
##    stock                             species        location               
##    <chr>                             <chr>          <chr>                  
##  1 Atlantic cod - Georges Bank       Atlantic cod   Georges Bank           
##  2 Atlantic cod - Gulf of Maine      Atlantic cod   Gulf of Maine          
##  3 Haddock - Georges Bank            Haddock        Georges Bank           
##  4 Haddock - Gulf of Maine           Haddock        Gulf of Maine          
##  5 Windowpane - Gulf of Maine / Geo… Windowpane     Gulf of Maine / George…
##  6 Windowpane - Southern New Englan… Windowpane     Southern New England /…
##  7 Winter flounder - Gulf of Maine   Winter flound… Gulf of Maine          
##  8 Winter flounder - Southern New E… Winter flound… Southern New England /…
##  9 Yellowtail flounder - Cape Cod /… Yellowtail fl… Cape Cod / Gulf of Mai…
## 10 Yellowtail flounder - Georges Ba… Yellowtail fl… Georges Bank           
## 11 Yellowtail flounder - Southern N… Yellowtail fl… Southern New England /…
## 12 Winter flounder - Georges Bank    Winter flound… Georges Bank           
## 13 Red hake - Gulf of Maine / North… Red hake       Gulf of Maine / Northe…
## 14 Red hake - Southern Georges Bank… Red hake       Southern Georges Bank …
## 15 Silver hake - Gulf of Maine / No… Silver hake    Gulf of Maine / Northe…
## 16 Silver hake - Southern Georges B… Silver hake    Southern Georges Bank …
## 17 Atlantic cod - Eastern Georges B… Atlantic cod   Eastern Georges Bank   
## 18 Haddock - Eastern Georges Bank    Haddock        Eastern Georges Bank   
## 19 Goosefish - Gulf of Maine / Nort… Goosefish      Gulf of Maine / Northe…
## 20 Goosefish - Southern Georges Ban… Goosefish      Southern Georges Bank …

Looks like there are just 8.

2 Gapfill

Use a linear model to fill years between assessments.

data_gf <- data %>%
  group_by(stock, indicator) %>%
  mutate(timeseries = n()) %>% #get number of years for each stock/indicator combo. We need at least 2 years
  filter(timeseries > 1) %>%
  complete(year = 2005:2018) %>%
  ungroup() %>%
  mutate(value2 = ifelse(year == 2018 & is.na(value), zoo::na.locf(value), value), #if 2018 is NA, use 2017 (or last reported value)
         value3 = ifelse(year == 2005 & is.na(value2), zoo::na.locf(value2, fromLast = TRUE), value2)) %>% #if 2005 is NA, roll back the closest value from following years
  group_by(stock, indicator) %>%
  mutate(gf_value = zoo::na.approx(value3), #linear interpolation to gapfill
         gf = ifelse(is.na(value), 1, 0)) %>% #if a value is gapfilled, assign it a 1, otherwise 0
  select(-value, -value2, -value3, -timeseries) %>%
  rename(value = gf_value)

#grab just B/Bmsy data
  nmfs_b_bmsy <- data_gf %>%
    filter(indicator == "b_bmsy") %>%
    select(stock, year, value)

#grab just F/Fmsy data  
  nmfs_f_fmsy <- data_gf %>%
    filter(indicator == "f_fmsy") %>%
    select(stock, year, value)
  ### 
  b_bmsy_underexploit_penalty <- 0.25
  b_bmsy_underexploit_thresh  <- 3.00
  f_fmsy_underfishing_penalty <- 0.25
  f_fmsy_overfishing_thresh   <- 2.00
  
  ### Apply rolling mean to F/Fmsy
  ## Why do we do this? because B is a less sensitive metric (relies of biological processes) and F can fluctuate pretty easily because it is really just a mgmt decision.
  nmfs_f_fmsy <- nmfs_f_fmsy %>%
    arrange(stock, year) %>%
    group_by(stock) %>%
    mutate(f_fmsy_rollmean = zoo::rollmean(value, k = 3, align = 'right', fill = NA)) %>%
    ungroup() %>%
    select(-value) %>%
    rename(value = f_fmsy_rollmean)
  
  stock_status_layers <- nmfs_b_bmsy %>%
    full_join(nmfs_f_fmsy) %>%
    spread(indicator, value) 
  
########################################################.
##### run each fishery through the Kobe plot calcs #####
########################################################.
### * ram_b_bmsy, ram_f_fmsy
  
  
### Function for converting B/Bmsy values into a 0 - 1 score
  rescale_bprime_crit <- function(fish_stat_df,
                                  bmax, bmax_val) {
    
    ###using NOAA's limits here
    overfished_th  <- 0.8
    ### 
    underfished_th <- 1.2
    
    bmax_adj <- (bmax - underfished_th) / (1 - bmax_val) + underfished_th
    ### this is used to create a 'virtual' B/Bmsy max where score drops
    ### to zero.  If bmax_val == 0, this is bmax; if bmax_val > 0, bmax_adj
    ### extends beyond bmax, to create a gradient where bmax_val occurs at bmax.
    
    fish_stat_df <- fish_stat_df %>%
      # group_by(stock) %>% ### grouping by stock will set b_max by max per stock, instead of max overall
      mutate(b_max     = max(b_bmsy, na.rm = TRUE)) %>%
      ungroup() %>%
      mutate(bPrime = NA,
             bPrime = ifelse(b_bmsy < overfished_th,  ### overfished stock
                             b_bmsy / overfished_th,
                             bPrime),
             bPrime = ifelse(b_bmsy >= overfished_th & b_bmsy < underfished_th,
                             1,                       ### appropriately fished stock
                             bPrime),
             bPrime = ifelse(b_bmsy >= underfished_th,
                             (bmax_adj - b_bmsy) / (bmax_adj - underfished_th), ### underfished stock
                             bPrime),
             bPrime = ifelse(bPrime < 0, 0, bPrime))
    return(fish_stat_df)
  }
  
  
  ### Function to create vertical gradient based on distance from
  ### ideal F/Fmsy value to actual F/Fmsy value
  f_gradient <- function(f, over_f, under_f, fmax, fmin_val) {
    x <- ifelse(f < over_f & f > under_f, 1, NA)
    x <- ifelse(f <= under_f, (f * (1 - fmin_val) / under_f + fmin_val), x)
    x <- ifelse(f >= over_f,  (fmax - f) / (fmax - over_f), x)
    x <- ifelse(f > fmax, NA, x)
    return(x)
  }
  
  ### Function to convert F/Fmsy values into 0 - 1 score
  rescale_fprime_crit <- function(fish_stat_df,
                                  fmax, fmin_val) {
    
    ### params - taken from BC but changed Bcrit to 0 instead of 0.4:
    Bcrit <- 0.5; overfished_th <- 0.8
    ### underfishing_th is set to the idea "1/3 for the birds":
    underfishing_th <- 0.66; overfishing_th  <- 1.2
    
    bcritslope = 1 / (overfished_th - Bcrit)
    ### connecting from (Bcrit, 0) to (overfished_th, 1)
    
    fish_stat_df <- fish_stat_df %>%
      mutate(fPrime = ifelse(b_bmsy < overfished_th & f_fmsy < fmax,
                             f_gradient(f_fmsy + (overfished_th - b_bmsy) * bcritslope,
                                        over_f = overfishing_th,
                                        under_f = underfishing_th,
                                        fmax = fmax,
                                        fmin_val = fmin_val),
                             NA),
             fPrime = ifelse(b_bmsy >= overfished_th & f_fmsy < fmax,
                             f_gradient(f_fmsy,
                                        over_f = overfishing_th,
                                        under_f = underfishing_th,
                                        fmax = fmax,
                                        fmin_val = fmin_val),
                             fPrime),
             fPrime = ifelse(is.na(fPrime), 0, fPrime), ### fill zeros everywhere unscored
             fPrime = ifelse(is.na(f_fmsy), NA, fPrime) ### but if no f_fmsy, reset to NA
      )
    return(fish_stat_df)
  }
  
  stock_status_df <- stock_status_layers %>%
    rescale_bprime_crit(bmax     = b_bmsy_underexploit_thresh,
                        bmax_val = b_bmsy_underexploit_penalty) %>%
    rescale_fprime_crit(fmax     = f_fmsy_overfishing_thresh,
                        fmin_val = f_fmsy_underfishing_penalty) %>%
    mutate(x_prod = ifelse(!is.na(fPrime), (fPrime * bPrime), bPrime),
           basis  = case_when(
             !is.na(fPrime) & !is.na(bPrime) ~ 'F_Fmsy, B_Bmsy',
             is.na(fPrime)  & !is.na(bPrime) ~ 'B_Bmsy only',
             is.na(bPrime)  & !is.na(fPrime) ~ 'F_Fmsy only'
           )) %>%
    dplyr::select(year, stock,
                  score = x_prod,
                  basis,
                  bPrime, fPrime,
                  b_bmsy, f_fmsy) 

stock_scores <- stock_status_df %>%
  filter(basis != "F_Fmsy only") %>%
  group_by(stock, year) %>%
  mutate(score = mean(score, na.rm = T)) %>%
  select(year, stock, score, bPrime, fPrime, b_bmsy, f_fmsy, basis) 

ggplot(stock_scores, aes(x = year, y = score)) +
  geom_line() +
  theme_bw() +
  facet_wrap(~stock)

kobe_plot <- function(sp){
  
 ss_df <- stock_scores %>%
    filter(stock == sp) %>%
    arrange(year) %>%
    mutate(last_bbmsy = last(b_bmsy),
           last_ffmsy = last(f_fmsy),
           last_datayear = last(year)) %>%
   ungroup()

generate_kobe_df <- function(f_fmsy_max = 2.5,
                             b_bmsy_max = 3.0,
                             reso       = 0.01,
                             bmax_val   = 0,
                             fmin_val   = 0,
                             weighting_b = 1) {

  kobe_raw <- data.frame(stock  = 1,
                     f_fmsy = rep(seq(0, f_fmsy_max, reso), each  = round(b_bmsy_max/reso) + 1),
                     b_bmsy = rep(seq(0, b_bmsy_max, reso), times = round(f_fmsy_max/reso) + 1))

  kobe <- kobe_raw %>%
    rescale_bprime_crit(bmax = b_bmsy_underexploit_thresh,
                        bmax_val = bmax_val) %>%
    rescale_fprime_crit(fmax = f_fmsy_overfishing_thresh,
                        fmin_val = fmin_val) %>%
    mutate(x_geom  = (fPrime * bPrime),
           x_arith = (fPrime + bPrime) / 2)

  return(kobe)
}

bbmsy_lim <- max(round(max(ss_df$b_bmsy, na.rm = TRUE) + .1, 1), 3)
ffmsy_lim <- max(round(max(ss_df$f_fmsy, na.rm = TRUE) + .1, 1), 2.5)
  
kobe_df <- generate_kobe_df(f_fmsy_max = ffmsy_lim,
                           b_bmsy_max = bbmsy_lim,
                           bmax_val = .25,
                           fmin_val = .25)


kobe_stock_plot <- ggplot(data = kobe_df, aes(x = b_bmsy, y = f_fmsy)) +
    theme_bw() +
    geom_raster(alpha = .8, aes(fill = x_geom)) +
    scale_fill_distiller(palette = 'RdYlGn', direction = 1) +
    labs(title = as.character(sp),
         x = 'B/Bmsy',
         y = 'F/Fmsy',
         fill = "Stock score") +
    geom_path(data = ss_df, 
              show.legend = FALSE,
              aes(x = b_bmsy, y = f_fmsy, group = sp),
              color = 'grey30') +
    geom_point(data = ss_df, 
               show.legend = FALSE,
              aes(x = last_bbmsy, y = last_ffmsy)) +
    geom_text(data = ss_df %>%
                mutate(year = ifelse(year/5 == round(year/5) | year == last_datayear, year, NA)), 
              aes(x = b_bmsy, y = f_fmsy, label = year), 
              hjust = 0, nudge_x = .05, size = 2)

return(kobe_stock_plot)
}

sp <- stock_scores %>%
  filter(basis == "F_Fmsy, B_Bmsy")

plots <- lapply(unique(sp$stock), kobe_plot)
cowplot::plot_grid(plotlist = plots[1:4])

cowplot::plot_grid(plotlist = plots[5:8])

cowplot::plot_grid(plotlist = plots[9:12])

cowplot::plot_grid(plotlist = plots[13:16])

cowplot::plot_grid(plotlist = plots[17:20])

cowplot::plot_grid(plotlist = plots[21:24])

cowplot::plot_grid(plotlist = plots[25:28])

cowplot::plot_grid(plotlist = plots[29:32])

cowplot::plot_grid(plotlist = plots[33:35])

write.csv(stock_scores, file = "data/nmfs_stock_scores.csv")